function y = D_inv_demand(i_m,D_grid,func_var)
%generate inverse demand function and the derivative of inverse demand for D_grid
D_grid = D_grid(:);
alpha1 = func_var.alpha1; alpha2 = func_var.alpha2; alpha3 = func_var.alpha3;
lambda_func = func_var.lambda_func; dlambda_func=func_var.dlambda_func;
NDgrid = length(D_grid);p_bar = func_var.p_bar;beta=func_var.beta;

z = zeros(NDgrid,1);
for i =1:length(D_grid)
    if alpha1*lambda_func(0)+alpha3*lambda_func(0+D_grid(i))-i_m<0
       z(i)=0;
    else
    if i<=5
        z(i) = fzero(@(z)alpha1*lambda_func(z)+alpha3*lambda_func(z+D_grid(i))-i_m,[0,p_bar]);
    else
        z(i) = fzero(@(z)alpha1*lambda_func(z)+alpha3*lambda_func(z+D_grid(i))-i_m,[0,z(i-1)+0.01]);
    end
  
    end
end
psi = beta*(alpha2*lambda_func(D_grid)+alpha3*lambda_func(z+D_grid))+beta;
A  = dlambda_func(z+D_grid);
B  = dlambda_func(z);
C = dlambda_func(D_grid);
dzdD=-alpha3*A./(alpha3*A+alpha1*B);
dpsi = beta*(alpha2*C+alpha3*A.*(1+dzdD));

y = [D_grid,z,psi,dpsi];